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Abstract 

Incompressible two-dimensional flows such as the advection (Li- 
ouville) equation and the Euler equations have a large family of 
conservation laws related to conservation of area. We present two 
Eulerian numerical methods which preserve a discrete analog of 
area. The first is a fully discrete model based on a rearrange- 
ment of cells; the second is more conventional, but still preserves 
the area within each contour of the vorticity field. Initial tests 
indicate that both methods suppress the formation of spurious 
oscillations in the field. 

1 Introduction 

When a smooth field u)(x,y) is advected by an area-preserving flow, the area 
within each contour of ui is preserved. This is seen in pure advection and in 
the Euler equations, for example, and is important in the numerical solution of 
two-phase free boundary problems, where the total volume of each fluid should 
be preserved. Yet, although the advection problem has been addressed in prob- 
ably thousands of papers, and very accurate, stable, and efficient methods are 
known, no existing numerical methods take the area-preservation property into 
account. In this Letter we present an initial study containing two methods 
which do preserve (a discrete analog of) area. Although they are not, presum- 
ably, competitive with the best existing methods for advection, the results are 
extremely promising. 

The configuration space of an inviscid incompressible fluid is the group T> M of 
volume-preserving diffeomorphisms of the fluid's domain; the 'Arnold' picture, 
in which the Euler equations are geodesic equations on this group equipped with 
the kinetic energy (L2) metric, is treated in p|. The configuration at any time 
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is a volume-preserving rearrangement of the initial condition. Existing Eulerian 
numerical methods do not preserve any discrete analogue of this property. This 
is particularly relevant in two dimensions, where area preservation leads to an 
infinite number of conserved quantities, the generalized enstrophies. 

We consider a two-dimensional fluid with divergence-free velocity field u = 
(u, v), stream function ip (i.e. u = ijj y , v = —ip x ), and some quantity w, which 
we call the vorticity, which is advected by the fluid: 

w + u- Vu; = ib + J[w, = 0, (1) 

where the Jacobian 

<*> 

In the two situations we shall consider, this is a Hamiltonian system with Poisson 
bracket 

{F, G} = J ujJ{6F, 6G) dxdy (3) 

and Hamiltonian 

1 f 

2i ^ w dxdy 

where the stream function ip is either a given function ip = ip(x,y,t), in which 
case Eq. (Q) is the Liouville (advection) equation, or is determined by the 
Poisson equation V 2 ip = —lu, in which case (|l|) is the the two-dimensional Euler 
equation. Other two-dimensional flows such as the shallow water and semi- 
geostrophic equations also possess a quantity u>, called a potential vorticity, 
that is advected according to Eq. There are also applications to level-set 
methods, in which to is not a physical variable but is introduced so that the 
curve u)(x,y) = c can indicate a free boundary. 

The Casimirs of the Poisson bracket (||) are conserved quantities of the PDE 
(§. These can be variously written as 



Cf = J f(u) dxdy 
for any function / such that Cf exists, as 

C„ = / uj u dxdy, 



called the generalized enstrophies (Ci is the usual enstrophy), or as the areas 
enclosed by each vorticity contour 

A(c) = 1 dxdy. 

J UJ>C 

They all reflect the fact that oj is being advected by an area-preserving vector 
field and can only reach states which are area-preserving rearrangements of its 
initial state. That is, 

uj(x,t) = ui(<pt 1 (x),Q), 
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where ipt is the time-i flow of the vector field u. 

The famous Arakawa Jacobian is an Eulerian finite difference approximation 
of Eq. (0) which preserves discrete analogues of the energy H and the enstrophy 
C 2 §7 ft is known to preserve the mean of the energy spectrum and to prevent 
some nonlinear instabilities. However, the other conserved quantities are not 
preserved and their role in the dynamics is not known [0. 

Area preservation can also be studied in a Lagrangian framework — for exam- 
ple, point vortex methods could be said to be area-preserving — but Lagrangian 
schemes carry a lot of extra information (the particle paths) which should be 
decoupled from the dynamics. The dimension of the phase space is halved in 
Eulerian form, and further reduced by preserving (discrete analogues of) the 
Casimirs. For ODEs, it is well established that the best long-time results are 
obtained by working in the smallest possible phase space ||. 

The Hamiltonian picture has been described by Marsden and Weinstein Q . 
The configuration space is the group V^. The Euler equations in Lagrangian 
form are a canonical Hamiltonian system on T*2? M , and in Eulerian form are 
a Lie-Poisson system on the dual of the Lie algebra of 2?^, which is identified 
with the space of vorticities. The coadjoint orbits of this space are the level 
sets of the Casimirs, each of which is a symplectic manifold. Discretizations 
of the Eulerian form are not, in general, Hamiltonian systems, nor do they 
have conserved quantities corresponding to the Casimirs (although there is one 
interesting Hamiltonian discretization, the sine- Euler equations JlO[ ) . 

Therefore we forget about the Hamiltonian structure and study the Casimirs — 
the area-preservation — and present two models in which a discrete analogue of 
the areas A(c) is preserved. The first (Section 2), based on a literal rearrange- 
ment of cells, is interesting in that it gives a fully-discrete, cellular-automata-like 
model of an incompressible fluid, ft does not preserve smoothness of the vor- 
ticity field (although filamentation and turbulence mean that it can't usually 
stay very smooth anyway). A smooth version (Section 3) is based on computing 
an approximation of A(c) which is smooth as a function of c, and relabelling 
the vorticity field so that A(c) is constant in time. It is tested on the Liouvillc 
equation and prevents the appearance of large spurious maxima and minima in 
the vorticity field during its evolution. 



2 The cell rearrangement model 

Both of the models presented here are projection schemes. The vorticity is 
evolved by any sensible scheme for some short time t (e.g., 1-10 time steps), 
and then projected onto some space of rearrangements of the original vorticity. 

In this section we consider the vorticity field to be piecewise constant on a 
set of fixed cells, which for convenience we take to be squares with side h. A 
(minuscule!) subset of the rearrangements of the initial condition is given by 
the permutations of the cells. However, these can be naturally associated with 
the fluid flow. For, consider the area-preserving map if which is the time-i flow 
of the fluid. According to a theorem of Lax H, there is a mapping P which 
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permutes cells and which satisfies 

1. P(C) n cp(C) ^ for all cells C; and 

2. \\P(x) - <p(x)\\ < sup VjZ&c My) - <p(z)\\ +V2h Vie C. 

The dynamics of such lattice maps are often studied. For example, if the con- 
tinuous map cp is iterated on a computer, it will not be exactly a bijection or 
exactly area-preserving, due to round-off error. By replacing it with a lattice 
map and examining the limit h — ► the effects of roundoff error can be studied. 

The easiest way to construct lattice maps is as a composition of shears 
x\ = Xi for i — 1, . . . ,d, x'j — \_fj(xi, . . . , Xd)\ for j = d + 1, . . . , n, where 
\x\ is the nearest lattice point to x. This would be suitable, for example, if ip 
itself were approximated by a product of shears, as is for example the flow of 
separable Hamiltonians H = H i (p) + H2 (q) (the flow of the Hamiltonian vector 
field corresponding to each Hi is a shear) . This is very fast and the permutation 
need not be constructed explicitly. 

However, in the present case ip can only be obtained by integrating the 
Lagrangian particle paths for a short time t, and an explicit approximating 
lattice map seems to be unobtainable. Scovel M suggested using maps of the 
form, e.g., x' = x+[J\7 S((x+x')/2)\ for a suitable Poincare generating function 
S (here S = (At)tp would give a good approximation of the time- At flow of the 
stream function tjj). However, this nonlinear, discrete equation does not seem 
to have solutions in general. 

Thus, it seems that one must laboriously construct a table of the permu- 
tation. An algorithm which does this is described in |Q. Its running time 
is 0(N 3 ), where N — 0{l/h 2 ) is the number of cells. One must construct 
lists of candidate cells (e.g., all those that intersect f(C)) and make successive 
choices from these lists, backtracking when no choices remain. While practical 
for moderate N when the dynamics of the lattice map are going to be studied 
intensively, in the present application ip changes at every time step; searching 
for a completely new permutation every time is too expensive. This approach 
has been explored by Turner Q . 

Luckily, there is a way out of this impasse, using the extra physical informa- 
tion attached to each cell: the vorticity itself. The only use of the permutation 
P is to update the vorticity field w by w h> w, 5 o F = w, in order that the 
distribution of vorticity values remains constant. This can be achieved directly, 
without actually constructing a P which approximates ip, by the following algo- 
rithm. Let rankj (c) be the number of cells with vorticities greater than or equal 
to lo at time t, i.e., 

rank t (c) = #{j : oj{xj,t) > c}, 

with ties broken arbitrarily to make rank t an invertible function onto {1, . . . , N}. 
Then: 

1. Update the field u> for time t any standard Eulerian method; and 

2. let ujj — rank^^rank^ojj)). 
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The new field uj can be constructed in time NlogN by sorting the two lists of 
vorticity values at times and t. The largest current value is replaced by the 
largest original value, and so on. (Other updates, based on minimizing \\ui — w||, 
are also possible.) 

This algorithm can be regarded as constructing a permutation, albeit a per- 
mutation that has no relationship to the flow ip. This is a truly finite-state model 
of an incompressible fluid: the state space is the permutation group Sn and the 
fluid dynamics reduces to the dynamics of the map Sn — * Sn, defined above, 
parameterized by the initial distribution of vorticity values. It is an almost 
cellular-automata-likc fluid model, although lacking the local update property 
of CAs. It has the aesthetic appeal of capturing the vorticity-rearrangement 
property perfectly in a naturally discrete way, of constructing a "discrete coad- 
joint orbit" , and it is very cheap. 

However, these advantages are offset by a practical disadvantage of lack of 
smoothness. The new vorticity values are selected somewhat arbitrarily from 
the sorted list, and the new field may be rougher than the original. This is 
probably unavoidable, given the chosen fully discrete state space. The noise 
of this imposed roughness may swamp any gains from preserving the coadjoint 
orbits. However, in a turbulent flow with highly filamented vorticity, the loss 
of smoothness may not be significant. A second consequence of the lack of 
smoothness is that if \u(x,t) — uj(x, 0)| is too small, then uj = u> — the field 
cannot be updated at all. The remapping interval t must be large enough to 
allow some change in the configuration. For example, the flow map ip should 
move each cell across at least 2 cells so that the algorithm has some scope for 
finding a suitable permutation. 

3 The vorticity relabelling model 

The cell rearrangement model produces an area function A(c) which is discon- 
tinuous — in fact, it is piecewise constant. To improve it, we need to 

(i) produce a smoother approximation of A(c). If ui(x, y) is a smooth function, 
we want an approximation of A(c) which is as smooth and accurate as 
possible, using only the grid values w(xi 7 yj); and 

(ii) project the vorticity function so that its area function A(c) at time t > 
equals (or closely approximates) the initial area function. 

3.1 Computing the areas enclosed by vorticity contours 

We consider a compact domain f2 with area 1, usually a square or torus, on 
which u) is bounded with range [w m i n , w max ], and of smoothness C r . It may be 
degenerate, e.g., constant on open sets. 

Definition 1 The area function of the field uj is the area enclosed by the set 
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{(x,y) : u(x,y) > c}, i.e., 



A u : [Wmin, W, 



max 




uj(x)>c 



Idxdy 



A(c) is strictly decreasing with respect to c. It is C r at regular (noncritical) 
values c, C° at nondegenerate critical values, and discontinuous at c if the set 
{x : oj(x,t) = c} has positive area. (Lack of differentiability at critical values 
can be seen by studying u = —{x 2 + y 2 ), for which A(c) = ire for c < and 
for c > 0.) Thus, its inverse A -1 exists and is nonincreasing (i.e., more area 
must be enclosed by a lesser value of uj.) 

Let X be an interpolation or approximation operator mapping grid functions 
to fields, i.e. functions on f2. 

Definition 2 The area function of the grid function ui is defined to be the area 
function of its interpolant, i.e., 



It automatically inherits the monotonicity properties of A. Let 1Z be a restriction 
operator mapping fields to grid functions, usually by evaluating on the grid. Let 
uj = ITZuj. The crucial observations are the following: 

1. Choice of 1 can lead to A~ being as smooth as A u , and of any order of 
accuracy as an approximation; 

2. If uj is piecewise linear, its contours are polygons, whose area can be found 
quickly for any contour topology; 

3. If uj has polygonal contours, A~ can be second order accurate and as 
smooth as A u . 

Item (1) is obvious, and is a consequence of existence of C r approximations to 
functions. An algorithm for finding the area enclosed by (unions of) polygons 
is given below. The most important point is (3), as it says that smooth inter- 
polants, which are expensive in two dimensions, are not needed to compute a 
smooth area function. 

Consider a grid function on a triangulation of Q. Interpolating by piecewise 
polynomials along edges only, and constructing the interpolant whose contours 
are line segments within each triangle (whose graph is a "ruled surface"), yields 
an area function which is as smooth as the interpolant at vertex (grid point) 
values and analytic elsewhere. Thus, only smooth one-dimensional interpolation 
is needed, which is relatively cheap (e.g., C 1 can be achieved using local cubics). 

Piecewise linear interpolation yields a C°, second-order-accurate area func- 
tion. However, it is better than its mere continuity might make it appear, since 
its derivative jumps at vertex values are only 0(h 2 ) on a grid with spacing h. 
So, numerically, it is indistinguishable from a C 1 function. In practice, the most 
glaring jumps are in its second derivative at vertex values, not in the function 
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itself. (See Fig. 2.) Piecewise linear interpolation seems to be suitable in prac- 
tice and this is what we use in the tests below. (If the main computational grid 
is square, we triangulate using an extra vertex at the center of each cell, whose 
value is assigned by linear interpolation.) However, true C 1 area functions have 
been tested as well. 

The great advantage of polygonal contours is that the area of a simple poly- 
gon with vertices Xj, i = 1, . . . , n, is very easy to compute: it is 



where x n+ \ := x\. This can be seen by deriving it for a triangle, triangulating 
the polygon against a fixed point, and then using independence with respect to 
the fixed point. It can also be viewed as a discretization of 



However, it would be expensive to chase contours around the domain and con- 
struct a list of simple polygons. Instead, one can simply scan each triangle for 
occurrence of a contour, find its endpoints xi, X2, and accumulate xi x X2 with 
a sign determined by the sense of the triangle when its vertices are visited in 
order of increasing function values. This handles arbitrary contour topology. 
(Exception handling is needed when two vertices and the contour all have equal 
values.) 

We are not sure if there is a similarly simple method with higher order con- 
tours. In practice, to get more than second order accuracy, we use Richardson 
extrapolation from a coarser grid. 

For a list of contour values, the above algorithm involves scanning the cells 
once and accumulating areas of the relevant contours. If iV c values of the area 
function are needed, and the grid size is O(h), then each cell will contain 0(hN c ) 
contours on average, so the computation takes time 0(N c /h). In practice, we 
take iV c = 0(l/h) and build a function table, which is later interpolated as 
needed. Thus computing the areas takes 0(l/h 2 ), i.e., it is linear in the number 
of grid points. 

If ui is nearly constant on large areas, then A{c) can be very steep, so it should 
be tabulated using adaptive stepping in c. An example of the C° estimate of 
A(c) given by piecewise linear interpolation is shown in Fig. 2, together with 
the piecewise constant estimate given by simply counting the number of vertices 
where w > c. Its (numerical) derivative indicates its smoothness. Some care 
must be taken when interpolating to maintain monotonicity. 

We have also computed smooth approximations of (c) for random uj fields 
whose contours have complicated topology. 

3.2 Projecting the the space of rearrangements 

After evolving u for a short time t with Eulerian method, we have two grid 
functions, the vorticity at time 0, w(0), and at time t, u{t). We wish to project 
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1 dx dy = — / d(xdy — ydx) = — 
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ui(t) so that it is an (approximation of) a rearrangement of w(0). The projection 
should be small and should not destroy smoothness. Traditional methods for en- 
forcing constraints, such as steepest descents, appear to be completely infeasible 
because of the global and sensitive dependence of A(c) on the vertex values of 
u). Our proposed method is a continuous version of the sorted- assignment used 
in the cell rearrangement model of Section 2. In words, we compute the area 
enclosed by the contour through each vertex value and replace it by the value 
that originally enclosed that much area. The contour shapes and topologies do 
not change: only the values associated with each contour change. 

Definition 3 The relabelling projection on grid functions cj, « uj(xi,t) is de- 
fined by u(t) i— > w(t), where 

A u (o)(u>i) =A u (t< j (uJi) (4) 

for each vertex i. 

It is well defined by monotonicity of A(c). It has an obvious continuum analog 
(replacing i by x in Eq. (0)), which if applied to every value of 10 taken by a 
smooth vorticity field, with oj(0) and u>(t) both C r , yields a new field uj that is 
C r away from critical points of loq and Ut and C° at such critical points. 

To compute a good approximation of this projection quickly, the current 
area function is tabulated and interpolated at the vertex values, and then A~, Q -> 
(which, of course, does not change during the run) is evaluated by interpolation. 
Of course, we do not have a true projection in that 5/5, because interpolation 
errors in the contours do change the contour shapes by a small amount when 
the vertex values are changed. We do not quite get A^( )( c ) = A^(t)( c ) f° r au 
c. However, these errors can be controlled independently of the discretization 
error in u>, for example, by using a higher order approximation of A. In a 
numerical test, one application of Richardson extrapolation to the areas enclosed 
by piecewise linear contours gave |A^(o) ~ ~ 10~ 4 on a relatively coarse 

20 x 20 grid. By contrast, without the relabelling projection, errors in the area 
function rapidly reach order 1. 

If the vorticity is evolved for a short time t, with a method of spatial order p, 
spatial errors dominate the error in the area function which are 0(th p ). Thus, 
with t = o(l), the projection only alters the field by o(h p ), and the overall 
method (after evolution and projection), is still consistent of the same order p. 
The projection cannot correct any errors in the shapes of the contours, but it 
can stop those errors growing further by propagation of the false distribution 
of vorticity values, which is particularly bad for the 2D Euler equations, where 
those values determine the velocity field itself. 

4 Numerical tests 

Here we illustrate some short tests to validate our approach and show that it is 
indeed possible to compute and preserve area in an Eulerian method. We use 
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a coarse (20 x 20) grid which barely resolves the solution, and a crude (second 
order) finite difference approximation to the spatial differences, in order to test 
whether the method can correct the large oscillations and area errors that result. 

We solve the Liouville equation in fl = [0, l] 2 with initial field u> = exp(— 45(x— 
|) 2 — 15 (y — \) 2 ) advected by the velocity field with stream function V = 
sin(7ra:) sin(7ry). (See Figure 1.) This velocity field has shear, so ui rapidly rolls 
up into a tight spiral, mimicking the filamentation of vorticity in the Euler 
equations. The spatial derivatives in Eq. (]]]) are approximated by the Arakawa 
Jacobian, which is second order and preserves discrete analogues of energy and 
enstrophy. Although the discrete enstrophy ^ cuf is preserved, this does not 
help the scheme preserve areas any better than (nonconservative) central differ- 
ences do. 

Particles at the maximum of lu have a period of about 0.75. We integrate 
with a second order method for 400 times steps of At = 0.003, or total time 1.2, 
during which this maximum rotates 1.6 times around the centre of the square 
fi. Spatial errors completely dominate the total error at t = 1.2. 

Without any projection, oscillations rapidly develop and the distribution of 
vorticity values is not maintained at all well (see Figure 3(a)). A large minimum 
of tp = —0.69 forms, next to a spurious local maximum of ip — 0.46. The initial 
maximum of 1 has not been preserved but has decayed to 0.87. The comparison 
between the initial and final area functions (see Figure 2) shows that the area 
within most vorticity contours is not preserved at all. 

The area-preserving methods both involve periodically remapping the vortic- 
ity. If this period is too short (e.g. one time step), then the cell rearrangement 
model cannot update the vorticity at all. If it is too long, then not only the 
area but also the topology of the level sets can alter, which can not be corrected 
by the present methods. Once a small island of vorticity has been created, for 
example, it must be advected by the flow. 

We first consider the cell rearrangement model of section 2. Suppose the 
remapping is applied every N r time steps. This needs a large N r to yield a 
reasonably smooth remapped vorticity field; but if N r is too large then spu- 
rious maxima can evolve which are not removed by the remapping. With no 
remapping, this maximum reaches 0.46. With N r = 50, it reaches 0.26. With 
N r = 20, there is no isolated spurious maximum, but oscillations start to ap- 
pear within the main island. These grow worse at N r = 10. Therefore, N r = 20 
seems a reasonable balance, and the final field is shown in Figure 3(b). In one 
remapping period, the central peak moves across about 2 cells. 

This remapping is very fast, but it does not maintain smoothness of to, as 
can be seen here. In fact, it is surprising that it works even as well as it does 
in this example. However, the lack of smoothness would not be a problem in 
problems involving poorly-resolved turbulent fields. 

We consider now the vorticity relabelling model of section 3. In this model 
we are free to decrease the remapping interval N r as desired: we still obtain 
smooth results with N r = 1, for example. As N r is decreased, the results 
progressively improve. For N r = oo, 20, 10, and 5, the peak of the spurious 
maximum is at u) — 0.46, 0.09, 0.02, and 0.002, respectively. (Because of its 
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smooth interpolation, it cannot completely eliminate this maximum, as the cell 
rearrangement model does.) Results for N r = 10 are shown in Figure 3(d). The 
final field is very smooth, considering the coarse 20 x 20 grid, and very plausibly 
represents an element of the original state composed with an area-preserving 
diffeomorphism. One contour of the exact solution (found by particle tracking) 
is shown in the background. The computed solution has clearly suffered far 
too much diffusion, a result of using diffusive, non-upwinded second differences 
to approximate the advection term. Nevertheless, it is impressive that such 
information can be extracted from the same method that produced Figure 3(a), 
by merely imposing some conservation laws. 

Finally, Figure 3(c), shows the vorticity relabelling model applied to an 
even simpler spatial discretization, namely ordinary central differences. It is in 
fact more accurate that the Arakawa Jacobian (Fig. 3(d)), being slightly less 
diffusive. Thus, preserving areas lets one use much simpler finite differences and 
still maintain smooth, non-oscillatory solutions. 

Any of the techniques presented here can be combined with a more sophis- 
ticated underlying Eulerian scheme. If we used a high-order, low-diffusion up- 
winding scheme, for example, then area errors would have been much less than 
in Fig. 3(a); but they would still increase over time. Applying the vorticity 
relabelling would still improve the solution. 

5 Discussion 

The methods discussed here take into account one large family of conservation 
laws. This possibility raises many questions. What is the effect of using these 
methods for very long times? What is their effect on other conservation laws 
such as energy and symplecticity? How well do they work on larger applica- 
tions such as the shallow water equations? (For level-set applications, a simpler 
update, adding a constant to the advected field so that the area inside one 
particular level set is preserved, may be preferable.) 

More theoretically, is it possible to regard the 'equal area' functions as defin- 
ing a discrete phase space in which consistent approximations can be directly 
derived, instead of using brute force modification of existing methods? While 
desirable, this looks difficult, since we are not projecting to any well-defined 
manifold. Consider the subset of +1 defined by 

^(uh) =A (c), i = l,...,N 2 . 

This does have dimension 1 in general, but is formidably curled up on itself. 
It may be better to think of the constrained phase space as the configurations 
lying within some small distance of a manifold of dimension R k ; , as we 
are enforcing one curve's worth of constraints. 
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Figure 1: Initial condition for the test problem in Section 4. The contours show 
level sets of w(0) = exp(— 45(x— |) 2 — 15(y— ^) 2 ). The arrows show the vector 
field corresponding to the stream function ip = sin(7nn) sin(7ry) by which lo is 
advected. 
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Figure 2: Numerical computation of the area enclosed within vorticity contours, (a): C° approximation 
to A w ( ) using piecewise linear interpolation (Section 3). Here ui(0) is the initial condition shown in Figure 
1. (b): Piecewise constant approximation to A^ro) by sorting the list of vorticity values (Section 2). (c): 
Area function A^(t) after evolving for time t = 1.2 with no area preservation with an enstrophy-preserving 
scheme (Section 4). In the vorticity relabelling projection, vorticity values are mapped from this curve 
back to (a), (d): Finite difference approximation to cL4^( )(c) / 'dc, showing that, although only C°, for 
numerical purposes it can be regarded as being differentiable. The kinks in this derivative are due to u> 
being set to zero on the boundary. 
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Figure 3: Results for the advection problem on Fig. 1 after 1.6 rotations about the center, (a): Arakawa 
differences with no area preservation. A large negative blob of vorticity forms and spawns a secondary 
positive blob. The dotted contour indicates the exact solution, (b): Arakawa differences with cell 
rearrangement applied every 20 time steps (At = 0.001). (c): Central differences with vorticity relabelling 
applied every 10 time steps, (d) Arakawa differences with vorticity relabelling applied every 10 time steps. 
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